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Abstract 

A brief introduction to spline functions and 5-splines, and specifically to monotone 
spline functions - with code in R and C and with some applications. 
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Note: This is a working paper which will be expanded/updated frequently. All suggestions 
for improvement are welcome. The directory gih.stat.ucla.edu/splines has a pdf version, the 
complete Rmd hie with all code chunks, the bib hie, and the R source code. 


1 Introduction 

To define spline functions we hrst define a finite sequence of knots T = {tj}, with t\ < ■ ■ ■ < t p , 
and an order m. In addition each knot tj has a multiplicity mj , the number of knots equal 
to tj. We suppose throughout that mj < m for all j. 

A function / is a spline function of order m for a knot sequence {tj} if 

1. / is a polynomial iTj of degree at most m — 1 on each half-open interval Ij = [tj,tj+ 1 ) 
for j = 1 ,■■■ ,p, 

2. the polynomial pieces are joined in such a way that T>_ f{tf) = D+ f(tj) for s = 
0,1, • • • , m — mj — 1 and j — 1, 2, • • • ,p. 

Here we use T)_ ; and T>\ for the left and right s -derivative operator. If nrij = m for some 
j, then requirement 2 is empty, if rrij = m — 1 then requirement 2 means 7 Tj(tj) = TTj + i(tj), 
i.e. we require continuity of / at tj. If 1 < r m,j < m — 1 then / must, be m — nij — 1 times 
differentiable, and thus continuously differentiable, at tj. 

In the case of simple knots (with multiplicity one) a spline function of order one is a step 
function which steps from one level to the next at each knot. A spline of order two is piecewise 
linear, with the pieces joined at the knots so that the spline function is continuous. Order 
three means a piecewise quadratic function which is continuously differentiable at the knots. 
And so on. 


2 












2 Basic Splines 


Alternatively, a spline function of order m can be defined as a linear combination of B-splines 
(or basic splines) of order m on the same knot sequence. A B-spline of order m is a spline 
function consisting of at most m non-zero polynomial pieces. A B-spline B jm is determined 
by the m + 1 knots tj<-< tj +m , is zero outside the interval [tj, tj +m ), and positive in the 
interior of that interval. Thus if tj = tj +m then Bj }in is identically zero. 

For an arbitrary finite knot sequence t\, • • ■ ,t p , there are p — m B-splines to of order m to 
be considered, although some may be identically zero. Each of the splines covers at most m 
consecutive intervals, and at most m — 1 different B-splines are non-zero at each point. 


2.1 Boundaries 

B-splines are most naturally and simply defined for doubly infinite sequences of knots, that 
go to Too in both directions. In that case we do not have to worry about boundary effects, 
and each subsequence of m + 1 knots defines a B-spline. For splines on finite sequences of p 
knots we have to decide what happens at the boundary points. 

There are B-splines for tj, ■ ■ ■ , tj +m for all j — 1, • • • ,p — m. This means that the first m — 1 
and the last m — 1 intervals have fewer than m splines defined on them. They are not part of 
what De Boor (2001), page 94, calls the basic interval. For doubly infinite sequences of knots 
there is not need to consider such a basic interval. 

If we had m additional knots on both sides of our knot sequence we would also have m 
additional B-splines for j — 1 — m, ■ ■ • ,0 and m additional B-splines for j — p — m +1, • • • ,p. 
By adding these additional knots we make sure each interval [tj, tj + ±) for j — 1, • • • ,p — 1 
has m B-splines associated with it. There is stil some ambiguity on what to do at t p , but we 
can decide to set the value of the spline there equal to the limit from the left, thus making 
the B-spline left-continuous there. 

In our software we will use the convention to define our splines on a closed interval [a, b} with 
r interior knots a < t\ < ■ ■ ■ < t r < b, where interior knot tj has multiplicity rrij. We extend 
this to a series of p — M + 2 m knots, with M = Y7j= i m j , by starting with m copies of a, 
appending m,j copies of tj for each j = 1, • • • , r, and finishing with m copies of b. Thus a 
and b are both knots with multiplicity m. This defines the extended partition (Schumaker 
(2007), p 116), which is just handled as any knot sequence would normally be. 

2.2 Normalization 

B-splines can be defined in various ways, using piecewise polynomials, divided differences, or 
recursion. The recursive definition, first used as a defining condition by De Boor and Hollig 
(1985), is the most convenient one for computational purposes, and that is the one we use. 
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But first, the conditions we have mentioned only determine the 5-spline up to a normalization. 
There are two popular ways of normalizing 5-splines. The iV-splines iVy m , a.k.a. the 
normalized B-splines j or order m, satisfies 


E^W = i- (i) 

i 

Note that in general this is not true for all t, but only for all t in the basic interval. 
Alternatively we can normalize to M-splincs, for which 


r+oo 


M jt m(t)dt = 


t'j+k 


Mj t m{t)dt = 1. 


There is the simple relationship 


( 2 ) 




tj+m 


m 




( 3 ) 


2.3 Recursion 

The recursive definition, due independently to Cox (1972) for simple knots and to De Boor 
(1972) in the general case, is 


or 


M j>m (t ) = 


L j-\-m 


-Mj ,„_,(()+ W 


tn 


t 


(j'Ji 


j+m 


( 4 ) 




t — I : 


m+j —1 




tj +m t 

tj+rn tj +1 


Nj+l,m— 1 (t) ■ 


( 5 ) 


A basic result in the theory of 5-splines is that the different 5-splines are linearly independent 
and form a basis for the linear space of spline functions (of a given order and knot sequence). 


3 Computations 

Before introducing our new C code we review some the approaches we have used in the past. 
This will also give us the opportunity to make some comparisons. 
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3.1 Low Order Splines 


The R code in lowSpline.R has three functions to compute splines of order one, two, and 
three. It does not acknowledge any boundary values, so only looks at R-splines on an interval 
spanned by interior knots. The formulas we use are 


N j, iO) 


1 if tj < x < tj+ 1 , 
0 otherwise 


AM*) 


x -tj 
tj+i—tj 
b'+2~ x 

tj+2-tj+l 


0 


if tj < x < tj. |_i, 
if t j+ 1 < x < t j+2 , ■ 
otherwise 


Nj,z(x) 


' ( x ~tj ) 2 

(b'+i — b')(b+2 — b) 

+ 2 ~x) (x-tj + 1 )(t j+3 -x) 

{tj+2~ b)(b+2 — b+i) (b+3 — b+i)(b+2 — b'+i) 

(b+3-j) 2 

(b+3-b+i)(b+3-b+2) 

0 


if tj < x < tj + 1 , 
if tj+i — ^ ^ / - 2 ■ 

if f i+2 < a; < tj +3 , 
otherwise 


In the example of Ramsay (1988) the knots are 0.0, 0.3, 0.5, 0.6, and 1.0 (in the example 0 
and 1 are endpoints of the interval, but we’ll just treat them as interior knots in a longer 
sequence). Also note that Ramsay computes M-splines, while we compute RAsplines. 

We start with the p — m = 5 — 2 = 3 R-splines of order 2. The basic interval is [0.3, 0.6]. 



x 
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Figure 1: Piecewise Linear Splines with Simple Knots 

There are only two 5-splincs of order 3 on these knots, and the basic interval is empty. 



X 

Figure 2: Piecewise Quadratic Splines with Simple Knots 

The programs also work for multiple knots. Consider the example from De Boor (2001), page 
92. The knots are 0, 1, 1, 3, 4, 6, 6, 6, and the order is three. The basic interval is [1,6]. 
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Figure 3: Piecewise Quadratic Splines with Multiple Knots 
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In De Boor (2001), p 89, we find another example in which there are only two distinct knots, 
an infinite sequence of zeroes, followed by an infinite sequence of ones. In this case there are 
only m different 5-splincs, restriction to [0,1] of the familiar polynomials 



for j — 0, ■ • • , m — 1. 
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Figure 4: Quadratic Bernstein Basis 


3.2 Using GSL 


The GNU Scientific Library (Galassi et al. (2016)) has 5-spline code. The function 


gslSplineO in R calls the compiled gslSpline. c, which is linked with the relevant code 
from libgsl. dylib. We use the Ramsay example again. The GSL implementation auto¬ 
matically adds the extra boundary knots for the extended partition, which makes the basic 
interval [0,1]. 

knots <- c(0,.3,.5,.6,l) 
order <- 3 

x<-seq(0 , 1 , length = 1001) 
h <- matrix (0, 1001, 6) 
for (i in 1 : 1001) 

h[i,] <- gslSpline (x[i], order, knots) 
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Figure 5: Piecewise Quadratic Splines using GSL 


3.3 Using Recursion 

We have previously published spline function code, using an R interface to C code, in De 
Leeuw (2015a). That code translated the Fortran code in an unpublished note by Sinha 
(????) to C. There are some limitations associated with this implementation. First, it is 
limited to extended partitions with simple inner knots. Second, the function to compute 
R-spline values recursively calls itself, using the basic relation (5), which is computationally 
not necessarily the best way to go. 

innerknots <- c( .3, .5, .6) 

degree <- 2 

lowknot <- 0 

highknot <- 1 

x<-seq(0 , 1 , length = 1001) 

h <- sinhaBasis (x, degree, innerknots, lowknot, highknot) 
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Figure 6: Piecewise Quadratic Splinesusing Recursion 


3.4 De Boor 

In this paper we implement the basic BSPLVB algorithm from De Boor (2001), page 111, for 
normalized D-splines. There are two auxilary routines, one to create the extended partition, 
and one that uses bisection to locate the knot interval in which a particular value is located 
(Schumaker (2007), p 191). The function bsplineBasis() takes an arbitrary knot sequence. 
It can be combined with extendPartitionO, which uses inner knots and boundary points 
to create the extended partion. 


3.5 Illustrations 

For our example, which is the same as the one from figure 1 in Ramsay (1988), we choose 
a = 0, b = 1, with simple interior knots 0.3, 0.5, 0.6. First the step functions, which have 
order 1. 

innerknots <- c(.3,.5,.6) 
multiplicities <- c(l,l,l) 
order <- 1 
lowend <- 0 
highend <- 1 

x <- seq (le-6, l-le-6, length = 1000) 

knots <- extendPartition (innerknots, multiplicities, order, lowend, highend)$knots 
h <- bsplineBasis (x, knots, order) 
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k <- ncol (h) 
par (mf row=c (2 , 2) ) 
for (j in 1 :k) { 

ylab <- paste ("B-spline" , formatC(j, digits 
plot (x, h[, j], type="l" , col = "RED", lwd 

} 


1, width = 2, format = "d")) 

3, ylab = ylab, ylim = c(0,l)) 
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Figure 7: Zero Degree Splines with Simple Knots 


Now the hat functions, which have order 2, again with simple knots. 

multiplicities <- c(l,l,l) 
order <- 2 

knots <- extendPartition (innerknots, multiplicities, order, lowend, highend)$knots 

h <- bsplineBasis (x, knots, order) 

k <- ncol (h) 

par (mfrow=c(2,3)) 

for (j in 1 :k) { 

ylab <- paste ("B-spline" , formatC(j, digits = 1, width = 2, format = "d")) 
plot (x, h[, j], type="l", col = "RED", lwd = 3, ylab = ylab, ylim = c(0,l)) 

> 
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X X 

Figure 8: Piecewise Linear Splines with Simple Knots 

Next piecewise quadratics, with simple knots, which implies continuous differentiability at 
the knots. This are the N-splines corresponding with the M-splines in figure 1 of Ramsay 
(1988). 

multiplicities <- c(l,l,l) 
order <- 3 

knots <- extendPartition (innerknots, multiplicities, order, lowend, highend)$knots 

h <- bsplineBasis (x, knots, order) 

k <- ncol (h) 

par (mfrow=c(2,3)) 

for (j in l:k) { 

ylab <- paste("B-spline" , formatC(j, digits = 1, width = 2, format = "d")) 
plot (x, h[, j], type="l", col = "RED", lwd = 3, ylab = ylab, ylim = c(0,l)) 

} 
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Figure 9: Piecewise Quadratic Splines with Simple Knots 


X 


If we change the multiplicities to 1, 2, 3, then we lose some of the smoothness. 

multiplicities <- c(l,2,3) 
order <- 3 

knots <- extendPartition (innerknots, multiplicities, order, lowend, highend)$knots 
h <- bsplineBasis (x, knots, order) 
k <- ncol (h) 
par (mfrow=c(3,3)) 
for (j in l:k) { 

ylab <- paste("B-spline" , formatC(j, digits = 1, width = 2, format = 
plot (x, h[, j], type="l", col = "RED", lwd = 3, ylab = ylab, ylim = 

} 


"d")) 

c(0,l)) 
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Figure 10: Piecewise Quadratic Splines with Multiple Knots 


4 Monotone Splines 

4.1 I-splines 

There are several ways to restrict splines to be monotone increasing. Since 5-splines are 
non-negative, the definite integral of a 5-spline of order m from the beginning of the interval 
to a value x in the interval is an increasing spline of order m + 1. Integrated 5-splines are 
known as I-splines (Ramsay (1988)). Non-negative linear combinations /-splines can be used 
as a basis for the convex cone of increasing splines. Note, however, that if we use an extended 
partition, then all /-splines start at value zero and end at value one, which means their 
convex combinations are the splines that are also probability distributions on the interval. 
To get a basis for the increasing splines we need to add the constant function to the /-splines 
and allow it to enter the linear combination with either sign. 
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4.1.1 Low Order I-splines 


Straightforward integration, and using (3), gives some explicit formulas. If we integrate the 
step functions we get the piecewise linear /-splines. 


0 

1 


if x < tj , 

if tj < x < tj + 1 , 

if x > tj + 1 - 


And if we integrate the piecewise linear 5-splines of order 1 we get piecewise quadratic 
/-splines. 



Mj 2 {t)dt 


0 

(x-tj) 2 

< (h+i -t i)( f j+2-q) 

2—tj _|_ {x-tj +1 )(tj +2 -x) 

tj+2—tj (tj + 2—tj)(tj+2—tj+l) 

1 


if x < tj, 
if tj < x < tj+i, 
if t j+ 1 < x < t j+2 , 
if x > t j+2 . 


Both sets of /-splines are plotted in the next two figures, using R functions from lowSpline. R. 



X 

Figure 11: Monotone Piecewise Linear Splines with Simple Knots 
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Figure 12: Monotone Piecewise Quadratic Splines with Simple Knots 


4.1.2 General Case 


Integrals of I-splines are most economically computed by using the formula first given by 
Gaffney (1976). If i is defined by tj + g-\ < x < tj + g then 


rx l 

Jxi 771 r _ 0 


It is somewhat simpler, however, to use lemma 2.1 of De Boor, Lyche, and Schumaker (1976). 
This says 


/ M jtTn (t)dt = N e,m+i(x) - N i,m+ l(a), 


J a 


t>j 


£>j 


If we specialize this to /-splines, we find , as in De Boor (1976), formula 4.11, 


j+r 

M j: m(t)dt = J2 N e,m+i(x) 

£=j 

for x < tj +r+ 1 - This shows that /-splines can be computed by using cumulative sums of 
/j-spline values. 

Note that using the integration definition does not give a natural way to define increasing 
splines of degree 1, i.e. increasing step functions. There is no such problem with the cumulative 
sum approach. 
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4.2 Increasing Coefficients 


As we know, a spline is a linear combination of 5-splines. The formula for the derivative 
of a spline, for example in De Boor (2001), p 116, shows that a spline is increasing if the 
coefficients of the linear combination of 5-splines are increasing. Thus we can fit an increasing 
spline by restricting the coefficients of the linear combination to be increasing, again using 
the 5-spline basis. 

It turns out this is in fact identical to using /-splines. If the 5-spline values at n points are 
in an n x r matrix H, then increasing coefficients (3 are of the form [3 = Sa + 7 e r , where S is 
lower-diagonal with all elements on and below the diagonal equal to one, where a > 0 , where 
e r has all elements equal to one, and where 7 can be of any sign. So H[3 = ( HS)a + 7 e n . 
The matrix Z = HS is easily and cheaply found in R by 1 — t(apply(h, 1, cumsum )). 

4.3 Increasing Values 

Finally, we can simply require that the n elements of H[3 are increasing. This is a less 
restrictive requirement, because it allows for the possibility that the spline is decresing 
between data values. It has the rather serious disadvantage, however, that it does its 
computations in n-dimensioiial space, and not in r-dimensional space, where r = M + m, 
which is usually much smaller than n. Software for the increasing-value restrictions has been 
written by De Leeuw (2015b). In this paper, however, we prefer the cumsumO approach. It 
is less general, but considerably more efficient. 

4.4 Illustrations 

We use the same Ramsay example as before, but now cumulatively. First we integrate step 
functions with simple knots, which have order 1, using isplineBasis (). The corresponding 
I-splines are piecewise linear with order 2 . 
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Figure 23: 

Now we integrate the hat functions, which have order 2, again with simple knots, to find 
piecewise quadratic I-splines of order 3. These are the functions in the example of Ramsay 
(1988). 
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Figure 13: Monotone Piecewise Linear Splines with Simple Knots 
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Finally, we change the multiplicities to 1, 2, 3, and compute the corresponding piecewise 
quadratic I-splines. 
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Figure 15: Monotone Piecewise Quadratic Splines with Multiple Knots 


5 Time Series Example 

Our first example smoothes a time series by fitting a spline. We use the number of births 
in New York from 1946 to 1959 (on an unknown scale), from Rob Hyndman’s time series 
archive at http://robjhyndman.com/tsdldata/data/nybirths.dat. 


5.1 B-splines 

First we fit R-splines of order 3. The basis matrix uses x equal to 1 : 168, with inner knots 
12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, and interval [1,168], 

innerknots <- 12 * 1:13 
multiplicities <- rep(l,13) 
lowend <- 1 
highend <- 168 
order <- 3 
x <- 1:168 

knots <- extendPartition (innerknots, multiplicities, order, lowend, highend)$knots 
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h <- bsplineBasis (x, knots, order) 
u <- lm.fit(h, births) 

res <- sum ((births - h%*°/oU$coefficients) ~2) 
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Figure 16: Monotone Piecewise Quadratic Splines with Simple Knots 
The residual sum of squares is 229.3835417745. 


5.2 I-splines 

We now fit the /-spline using the 5-spline basis. Compute Z = HS using cumsumO, and 
then y and Z by centering (substracting the column means). The formula is 

min SSQ (y — Za — ye n ) = min SSQ (y — Za). 

a>0,7 a>0 

We use pnnlsO from Wang, Lawson, and Hanson (2015). 

knots <- extendPartition (innerknots, multiplicities, order, lowend, highend)$knots 
h <- isplineBasis (x, knots, order) 
g <- cbind (1, h[,-1]) 
u <- pnnls (g, births, l)$x 
v <- g°/„*°/„u 
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Figure 17: Monotone Piecewise Linear Splines with Simple Knots 
The residual sum of squares is 288.4054982424. 


5.3 B-Splines with monotone weights 


Just to make sure, we also solve the problem 


min 

Pl<P2<-<Pp 


SSQ (y-Xfi), 


which should give the same solution, and the same loss function value, because it is just 
another way to fit /-splines. We use the lsi() function from Wang, Lawson, and Hanson 
(2015). 

knots <- extendPartition (innerknots, multiplicities, order, lowend, highend)$knots 

h <- bsplineBasis (x, knots, order) 

nb <- ncol (h) 

d <- matrix (0,nb-l ,nb) 

diag(d)=-l 

dfouter(1 :(nb-1), 1 :nb,function(i,j) (j - i) == 1) ]<— 1 
u<-lsi (h,births ,e=d,f=rep(0,nb-l) ) 
v <- h 1*1 u 
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Figure 18: Monotone Piecewise Quadratic Splines with Simple Knots 
The residual sum of squares is 288.4054982424, indeed the same as before. 


5.4 B-Splines with monotone values 


Finally we solve 


min SSQ (y — X(3) 

xi/3<-<</3 


using mregnnMO from De Leeuw (2015a), which solves the dual problem using nnls() from 
Wang, Lawson, and Hanson (2015). 

knots <- extendPartition (innerknots, multiplicities, order, lowend, highend)$knots 
h <- bsplineBasis (x, knots, order) 
u <- mregnnM(h, births) 
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Figure 19: Monotone Piecewise Quadratic Splines with Multiple Knots 

The residual sum of squares is 288.3210359866, which is indeed smaller than the /-splines 
value, although only very slightly so. 


6 Regression Example 

We also analyze a regression example, using the Neumann data that have been analyzed 
previously in Gih (1990), pages 370-376, and in De Leeuw and Mair (2009),pages 16-17. The 
predictors are temperature and pressure, the outcome variable is density. We use a step 
function monotone spline for temperature and a piecewise quadratic monotone spline for 
pressure. The lest squares problem is to fit 

SSQ (y - (qe n + iZiaq + H 2 a 2 )), 

where H i and H 2 are /-spline bases and aq and a 2 are non-negative. We do not transform 
the outcome variable density, to keep things relatively simple. 

data(neumann, package = "homals") 

knotsi <- c(0, 20, 40, 60, 80, 100, 120, 140) 

orderl <- 1 

knots2 <- c(0, 0,0,100,200,300,400,500,600,600,600) 
order2 <- 3 

hi <- isplineBasis(200-neumann[,1] , knotsl, orderl) 
h2 <- isplineBasis(neumann[,2] , knots2, order2) 
g <- cbind (1, hi, h2) 
u <- pnnls (g, neumann[,3], l)$x 
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ul <- u[l+l:ncol(hl)] 

u2 <- u[l+ncol(hl)+l:ncol(h2)] 

The next two plots give the transformations of the predictors. 



temperature 

Figure 20: Regression Example: Piecewise Constant Temperature 



pressure 

Figure 21: Regression Example: Piecewise Quadratic Pressure 

And finally, we give the residual plot, observed density versus predicted density. 
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7 Appendix: Code 

7.1 R code 

7.1.1 lowSpline.R 

ZbSpline <- function(x, knots, k = 1) { 

ZbSplineSingle <- function(x, knots, k = 1) { 
kO <- knots[k] 
kl <- knots[k + 1] 
if ((x > kO) kk (x <= kl)) { 
return (1) 

> 

return(O) 

> 

return(sapply(x, function(z) 

ZbSplineSingle (z, knots, k))) 

} 

LbSpline <- function(x, knots, k = 1) { 

LbSplineSingle <- function(x, knots, k = 1) { 
kO <- knots[k] 
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kl <- knots[k + 1] 
k2 <- knots[k + 2] 
fl <- function(x) 

(x - kO) / (kl - kO) 
f2 <- function(x) 

(k2 - x) / (k2 - kl) 
if ((x > kO) kk (x <= kl)) { 
return(f1 (x)) 

> 

if ((x > kl) kk (x <= k2)) { 
return(f 2(x) ) 

> 

return(O) 

} 

return(sapply(x, function(z) 

LbSplineSingle(z, knots, k))) 

} 

QbSpline <- function(x, knots, k = 1) { 

QbSplineSingle <- function(x, knots, k = 1) { 
kO <- knots[k] 
kl <- knots[k + 1] 
k2 <- knots[k + 2] 
k3 <- knots[k + 3] 
fl <- function(x) 

((x - kO) ~ 2) / ((k2 - kO) * (kl - kO)) 
f2 <- function(x) { 

terml <- ((x - kO) * (k2 - x)) / ((k2 - kO) * (k2 - kl)) 
term2 <- ((x - kl) * (k3 - x)) / ((k3 - kl) * (k2 - kl)) 


return (terml + term2) 

> 

f3 <- function(x) 



( (k3 

- x) 

~ 2) 

/ 

((k3 - 

kl) * (k3 - k2)) 

if 

((x 

> kO) 

kk 

(x 

<= 

kl)) 

{ 


return(f 1 (x) ) 





> 








if 

((x 

> kl) 

kk 

(x 

<= 

k2) ) 

{ 


return(f 2(x) ) 





> 








if 

((x 

/"“S 

CN 

A 

kk 

(x 

<= 

k3)) 

{ 


return(f 3(x)) 






> 

return(O) 
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return(sapply(x, function(z) 

QbSplineSingle(z, knots, k))) 

} 

IZbSpline <- function(x, knots, k = 1) { 

IZbSplineSingle <- function(x, knots, k = 1) { 
kO <- knots[k] 
kl <- knots[k + 1] 
if (x <= kO) { 
return (0) 

> 

if ((x > kO) kk (x <= kl)) { 
return ((x - kO) / (kl - kO)) 

> 

if (x > kl) { 
return (1) 

> 

} 

return(sapply(x, function(z) 

IZbSplineSingle (z, knots, k))) 

} 

ILbSpline <- function(x, knots, k = 1) { 

ILbSplineSingle <- function(x, knots, k = 1) { 
kO <- knots[k] 
kl <- knots[k + 1] 
k2 <- knots[k + 2] 
fl <- function(x) 

((x - kO) ~ 2) / ((kl - kO) * (k2 - kO)) 
f 2 <- 

function(x) 

((x - kO) / (k2 - kO)) + ((x - kl) * (k2 - x)) / ((k2 - kO) * (k2 - kl)) 
if (x <= kO) { 
return (0) 

> 

if ((x > kO) kk (x <= kl)) { 
return(f1 (x)) 

> 

if ((x > kl) kk (x <= k2)) { 
return(f 2(x) ) 

> 

if (x > k2) { 
return (1) 

> 
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return(O) 


} 

return(sapply(x, function(z) 
ILbSplineSingle (z, knots, k))) 


7.1.2 gslSpline.R 

dyn. loadC'gslSpline .so") 

gslSpline <- function (x, k, br) { 
nbr <- length (br) 
nrs <- k + nbr - 2 
res <- 
• C( 

"BSPLINE" , 
as.double (x), 
as.integer (k), 
as.integer (nbr), 
as.double (br), 

results = as.double (rep(0.0, nrs)) 

) 

return (res$results) 


7.1.3 sinhaSpline.R 

dyn.load ("sinhaSpline .so") 
sinhaBasis <- 

function (x, degree, innerknots, lowknot = min(x,innerknots), highknot = max(x, innerkr 
innerknots <- unique (sort (innerknots)) 
knots <- 

c(rep (lowknot, degree + 1), innerknots, rep(highknot, degree + 1)) 
n <- length (x) 

m <- length (innerknots) + 2 * (degree + 1) 
nf <- length (innerknots) + degree + 1 
basis <- rep (0, n * nf) 
res <- .C( 

"sinhaBasis", d = as.integer (degree), 

n = as.integer (n), m = as.integer (m), x = as.double (x), knots = as.double (knot 
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} 


) 

basis <- matrix (res$basis, n, nf) 
basis <- basis[ ,which(colSums (basis) > 0)] 
return (basis) 


7.1.4 deboor.R 

dyn.load ("deboor.so") 

checklncreasing <- function (innerknots, lowend, highend) { 
h <- .C( 

"checklncreasing" , 

as.double (innerknots), 

as.double (lowend), 

as.double (highend), 

as.integer (length (innerknots)), 

fail = as.integer (0) 

) 

return (h$fail) 

> 

extendPartition <- 
function (innerknots, 

multiplicities, 
order, 
lowend, 
highend) { 

ninner <- length (innerknots) 
kk <- sum(multiplicities) 
nextended <- kk + 2 * order 
if (max (multiplicities) > order) 
stop ("multiplicities too high") 
if (min (multiplicities) < 1) 
stop ("multiplicities too low") 
if (checklncreasing (innerknots, lowend, highend)) 
stop ("knot sequence not increasing") 
h <- 
.C( 

"extendPartition" , 
as.double (innerknots), 
as.integer (multiplicities), 
as.integer (order). 


as.integer (ninner), 
as.double (lowend), 
as.double (highend), 

knots = as.double (rep (0, nextended)) 

) 

return (h) 

} 

bisect <- 
function (x, 

knots, 

lowindex = 1 , 

highindex = length (knots)) { 

h <- .C( 

"bisect" , 

as.double (x), 

as.double (knots), 

as.integer (lowindex), 

as.integer (highindex), 

index = as.integer (0) 

) 

return (h$index) 

} 

bsplines <- function (x, knots, order) { 

if ((x > knots [length(knots)] ) || (x < knots[l])) 
stop ("argument out of range") 
h <- .C( 

"bsplines" , 

as.double (x), 

as.double (knots), 

as.integer (order), 

as.integer (length (knots)), 

index = as.integer (0), 

q = as.double (rep(0, order)) 

) 

return (list (q = h$q, index = h$ind)) 


bsplineBasis <- function (x, knots, order) { 
n <- length (x) 
k <- length (knots) 
m <- k - order 
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result <- rep (0, n * m) 

h <- .C( 

"bsplineBasis" , 

as.double (x) , 

as.double (knots), 

as.integer (order), 

as.integer (k), 

as.integer (n), 

result = as.double (result) 

) 

return (matrix (h$result, n, m)) 

> 

isplineBasis <- function (x, knots, order) { 
n <- length (x) 
k <- length (knots) 
m <- k - order 
result <- rep (0, n * m) 
h <- .C( 

"isplineBasis" , 

as.double (x), 

as.double (knots), 

as.integer (order), 

as.integer (k), 

as.integer (n), 

result = as.double (result) 

) 

return (matrix (h$result, n, m)) 


7.2 C code 

7.2.1 gslSpline.c 

#include <gsl/gsl_bspline.h> 

void BSPLINE(double *x, int *order, int *nbreak, double *brpnts, 
double ^results) { 

gsl_bspline_workspace *my_workspace = 

gsl_bspline_alloc((size_t)border, (size_t)*nbreak); 
size_t ncoefs = gsl_bspline_ncoeffs(my_workspace); 
gsl_vector ^values = gsl_vector_calloc(ncoefs); 
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gsl_vector ^breaks = gsl_vector_calloc((size_t)*nbreak); 
for (int i = 0; i < *nbreak; i++) 

gsl_vector_set(breaks, (size_t)i, brpnts[i]); 

(void) gsl_bspline_knots(breaks, my_workspace); 

(void)gsl_bspline_eval(*x, values, my_workspace); 

for (int i = 0; i < ncoefs; i++) results [i] = (values->data)[i]; 

gsl_bspline_free(my_workspace); 

gsl_vector_free(values); 

gsl_vector_free(breaks); 

} 

7.2.2 sinhaSpline.c 

#include <stddef.h> 

#include <stdio.h> 

#include <stdlib.h> 

double bs(int nknots, int nspline, int degree, double x, double *knots); 

int mindex(int i, int j, int nrow); 

void sinhaBasis (int *d, int *n, int *m, double *x, double *knots, 
double *basis) { 
int mm = *m, dd = *d, nn = *n; 
int k=mm-dd- 1, i, j, ir, jr; 
for (i = 0; i < nn; i++) { 
ir = i + 1; 

if (x[i] == knots [mm - 1]) { 

basis[mindex(ir, k, nn) - 1] = 1.0; 
for (j = 0; j < (k - 1); j++) { 
jr = j + 1; 

basis[mindex(ir, jr, nn) - 1] = 0.0; 

} 

} else { 

for (j = 0; j < k; j++) { 
jr = j + 1; 

basis[mindex(ir, jr, nn) - 1] = bs(mm, jr, dd + 1, x[i], knots); 

} 

} 

> 

> 

int mindex(int i, int j, int nrow) { return (j - 1) * nrow + i; } 
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double bs(int nknots, int nspline, int updegree, double x, double *knots) { 
double y, yl, y2, tempi, temp2; 
if (updegree == 1) { 

if ((x >= knots[nspline - 1]) && (x < knots[nspline])) 

Y = 1.0; 

else 

y = 0.0; 

} else { 

tempi = 0.0; 

if ((knots[nspline + updegree - 2] - knots[nspline -1]) >0) 
tempi = (x - knots[nspline - 1]) / 

(knots[nspline + updegree - 2] - knots[nspline - 1]); 

temp2 = 0.0; 

if ((knots[nspline + updegree - 1] - knots[nspline]) > 0) 
temp2 = (knots[nspline + updegree - 1] - x) / 

(knots[nspline + updegree - 1] - knots[nspline]); 
yl = bs(nknots, nspline, updegree - 1, x, knots); 
y2 = bs(nknots, nspline + 1, updegree - 1, x, knots); 
y = tempi * yl + temp2 * y2; 

> 

return y; 

> 

7.2.3 deboor.c 


#include <math.h> 

#include <stdbool.h> 

#include <stdlib.h> 

inline int VINDEX(const int) ; 

inline int MINDEX(const int, const int, const int); 

void checkIncreasing(const double *, const double *, const double *, 

const int *, bool *); 

void extendPartition(const double *, const int *, const int *, const int *, 

const double *, const double *, double *); 
void bisect(const double *, const double *, const int *, const int *, int *); 
void bsplines (const double *, const double *, const int *, const int *, int *, 
double *); 

void bsplineBasis (const double *, const double *, const int *, const int *, 

const int *, double *); 

void isplineBasis (const double *, const double *, const int *, const int *, 
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const int *, double *); 


inline int VINDEX(const int i) { return i - 1; } 

inline int MINDEX(const int i, const int j, const int n) { 
return (i - 1) + (j - 1) * n; 

} 

inline int IMIN(const int a, const int b) { 
if (a > b) return b; 
return a; 

} 

inline int IMAX( const int a, const int b) { 
if (a < b) return b; 
return a; 

} 

void checkIncreasing(const double *innerknots, const double *lowend, 

const double *highend, const int *ninner, bool *fail) { 

*fail = false; 

if Olowend >= innerknots[VINDEX(l)]) { 

*fail = true; 
return; 

> 

if (*highend <= innerknots[VINDEX(*ninner)]) { 

*fail = true; 
return; 

> 

for (int i = 1; i < *ninner; i++) { 

if (innerknots[i] <= innerknots[i - 1] ) { 

*fail = true; 
return; 

} 

> 

} 

void extendPartition(const double ^innerknots, const int ^multiplicities, 

const int border, const int *ninner, const double *lowend, 
const double *highend, double ^extended) { 

int k = 1 ; 

for (int i = 1; i <= *order; i++) { 
extended[VINDEX(k)] = *iowend; 
k++; 
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> 

for (int j = 1; j <= *ninner; j ++ ) 

for (int i = 1; i <= multiplicities[VINDEX(j)]; i++) { 
extended[VINDEX(k)] = innerknots[VINDEX(j)]; 
k++; 

} 

for (int i = 1; i <= *order; i++) { 
extended[VINDEX(k)] = *highend; 
k++; 

> 

} 

void bisect(const double *x, const double *knots, const int *lowindex, 
const int *highindex, int *index) { 
int 1 = *lowindex, u = *highindex, mid = 0; 
while ((u - 1) > 1) { 

mid = (int)floor((u +1) / 2) ; 
if (*x < knots[VINDEX(mid)]) 
u = mid; 

else 

1 = mid; 

> 

*index = 1; 
return; 

} 

void bsplines (const double *x, const double *knots, const int *order, 
const int *nknots, int *index, double *q) { 
int lowindex = 1, highindex = *nknots, m = border, j, jpl; 
double drr, dll, saved, term; 

double *dr = (double *)calloc((size_t)m, sizeof (double) ); 
double *dl = (double *)calloc((size_t)m, sizeof (double) ); 

(void) bisect(x, knots, &lowindex, &highindex, index); 
int 1 = *index; 
for (j = 1; j <= m; j++) { 
q[VINDEX(j)] = 0.0; 

> 

if (*x == knots[VINDEX(*nknots)]) { 
q[VINDEX(m)] =1.0; 
return; 

> 

q[VINDEX(1)] = 1.0; 

j = i; 

if (j >= m) return; 
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while (j < m) { 

dr[VINDEX(j)] = knots[VINDEX(1 + j)] - * x; 

dl[VINDEX(j)] = *x - knots[VINDEX(1 + 1 - j)]; 

jpi = j + i; 

saved = 0.0; 

for (int r = 1; r <= j; r++) { 
drr = dr[VINDEX(r)]; 
dll = dl[VINDEX(jpl - r)]; 
term = q[VINDEX(r)] / (drr + dll); 
q[VINDEX(r)] = saved + drr * term; 
saved = dll * term; 

} 

q[VINDEX(jpl)] = saved; 

j = jpi; 

> 

free(dr); 
free(dl); 
return; 


void bsplineBasis (const double *x, const double *knots, const int *order, 

const int *nknots, const int *nvalues, double *result) { 
int m = *order, 1=0; 

double *q = (double *)calloc((size_t)m + 1, sizeof(double) ); 
for (int i = 1; i <= *nvalues; i++) { 

(void)bsplines(x + VINDEX(i), knots, order, nknots, &1, q); 
for (int j = 1; j <= m; j++) { 

int r = IMIN(1 - m + j, *nknots - m); 
result[MINDEX(i, r, *nvalues)] = q[VINDEX(j)]; 

} 

> 

free(q); 
return; 

} 

void isplineBasis(const double *x, const double *knots, const int border, 

const int *nknots, const int *nvalues, double ^result) { 
int m = *nknots _ *order, n = *nvalues; 

(void)bsplineBasis(x, knots, order, nknots, nvalues, result); 
for (int i = 1; i <= n; i++) { 

for (int j = m - 1; j >= 1; j—) { 

result[MINDEX(i, j, n)] += result[MINDEX(i, j + 1, n)]; 

} 

> 
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return; 

} 
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